function y = func_2_nonmod_analytical(q, zeta, m, CONSTS)

    k0 = CONSTS.k0;
    a  = CONSTS.a;
    c = CONSTS.c;
    eps_a = CONSTS.eps_a;
    eps = CONSTS.eps;
    eta = CONSTS.eta;
    Z0 = 4*pi/c;
    
    Q = k0*a*q;
    
    besseli_m = besseli(m, Q, 1);
    besseli_mp1 = besseli(m+1, Q, 1);
    
    Um = besseli_mp1.^2 + ((eps_a^2)/(eps*abs(eta)))*besseli_m.^2;
      
    %y = ((1i*Z0*m^2*eps_a)/(pi*k0*a^2*(eps_a^2+abs(eps*eta)))) * ...
    %    ((exp(-k0*q.*abs(zeta)))./(sqrt(q.^2-eps_a)));
    
    y = 1i*Z0*((m^2)/(pi*a))*(eps_a/(eps*abs(eta)))*((besseli_m.^2)./(Q.*Um)) .* ... 
        exp(-k0*sqrt(eps/abs(eta))*q.*abs(zeta));

end